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1.  SUMMARY 


A  Thermal  Protection  System  (TPS)  is  responsible  for  protecting  a  spacecraft’s 
components  from  melting  due  to  high  re-entry  temperatures.  In  the  design  of  TPS,  both 
maximum  thermal  stress  and  minimum  natural  frequency  must  be  considered  due  to  the 
combined  thermo-acoustic  environment  inherent  in  high-speed  vehicle  applications.  In  this  work, 
a  multi-objective  structural  optimization  method  for  the  three-dimensional  acreage  TPS  design  is 
developed  using  an  Evolutionary  Structural  Optimization  (ESO)  algorithm.  The  static  control 
parameter  used  to  find  the  optimum  in  minimum  thermal  stress  design  is  modified  to  address  an 
irregular  mode-switching  phenomenon,  as  well  as  for  improving  the  modal  stiffness  in  dynamic 
analysis.  Two  objectives  are  optimized  simultaneously;  namely,  the  maximization  of  fundamental 
natural  frequency  and  the  minimization  of  maximum  thermal  stress.  The  proposed  modified 
control  parameter  is  demonstrated  on  the  design  of  a  metallic  TPS  using  the  method  of  weighted 
objectives.  The  results  are  then  compared  with  the  conventional  ESO  sensitivity  approach.  This 
work  concludes  by  applying  the  methodology  which  makes  use  of  both  topology  and  shape 
optimization  in  the  design  of  an  acreage  TPS. 
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2.  INTRODUCTION 


The  primary  challenges  that  must  be  addressed  to  enable  lower-cost  access  to  space  are 
weight,  reusability,  and  ease  of  maintenance.  The  current  methodology  for  the  design  of  space 
and  high  speed  air  vehicles  is  to  construct  a  low-temperature  load  bearing  structure  and  then  to 
apply  a  thermal  protective  layer  [also  known  as  a  Thermal  Protection  System  (TPS)]  whose 
purpose  is  to  shield  the  low  temperature  structure.  A  successful  TPS  will  not  only  perform  its 
primary  function  of  maintaining  the  underlying  vehicle  structure  within  acceptable  temperature 
limits,  but  must  also  be  durable,  operable,  cost  elfective,  and  lightweight.  By  its  very  nature,  TPS 
is  parasitic.  Current  research  efforts  in  “hot”  structures  are  investigating  the  feasibility  of  using 
advanced  structural  materials  that  can  provide  both  structural  rigidity  and  survivability  in  the 
high-temperature  environment.  No  conclusion  has  yet  been  reached  which  of  the  two  approaches 
will  prove  to  be  the  most  efficient.  Hence,  continued  research  into  parasitic  TPS  is  extremely 
relevant. 

Existing  TPS  consists  of  different  types  of  materials  that  are  distributed  over  the 
spacecraft.  These  various  thermal  insulators  consist  of  felt  blankets,  ceramic  blankets,  ceramic 
tiles,  carbon-carbon  leading  edges,  as  well  as  metallic  TPS.  Metallic  TPS  has  obvious  advantages 
and  disadvantages;  the  nature  and  response  of  metals  is  well  understood.  Furthermore,  the 
durability  and  fracture  toughness  of  metals  are  generally  much  greater  than  that  of  ceramics.  One 
disadvantage  of  metallic  systems  is  weight.  This  work  attempts  to  mitigate  this  effect  through 
weight  minimization. 

The  metallic  TPS  design  for  a  spacecraft  operating  in  extreme  environments  of  thermal 
and  acoustic  loading  is  of  significant  importance.1  Design  considerations  for  thermal  and 
acoustic  loading  conditions  tend  to  be  at  odds  with  one  another.  Floating  or  compliant  designs 
will  tend  to  reduce  thermal  stresses,  whereas  less  compliant  designs  increase  natural  frequency 
and  avoid  excitation  from  wide-band  random  engine  noise.  To  address  the  TPS  design  efficiently, 
optimization  methods  are  essential  both  from  the  standpoint  of  identifying  a  potentially  small 
feasible  region  and  driving  toward  low  weight  designs. 

Traditional  structural  topology  optimization  methods,  such  as  the  density-based 
method,2  the  Homogenization  method,3’  4  and  the  Evolutionary  Structural  Optimization  (ESO) 
method,5'  6  are  applicable  to  multi-objective  optimization  problems  concerned  with  static  or 
dynamic  problems.  Multi-objective  optimization  is  the  method  of  forming  a  solution  that 
satisfies  a  number  of  conflicting  objectives  in  the  “best”  way  possible.  That  is,  multi-objective 
optimization  is  the  generation  of  a  design  that  achieves  the  optimum  perfonnance  of  the  structure 
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while  consideration  is  given  to  various  criteria.  ’  The  solution  to  the  multi-objective  problem  is 
known  as  a  Pareto  optimum.8’9  Pareto  optima  are  not  unique  solutions,  but  rather  they  consist  of 
a  series  of  solutions.10 

This  paper  incorporates  the  multi-objective  optimization  problem  into  the  ESO 
algorithm  using  the  weighted  objectives  technique.  The  ESO  method  developed  by  Xie  and 
Steven  in  1992  slowly  removes  inefficient  material  from  a  structure,  and  the  residual  shape  of  the 
structure  evolves  as  an  optimum.5  The  ESO  procedure  does  not  address  the  mode-switching 
phenomenon  that  is  common  in  three-dimensional  dynamic  problems.  Mode-switching  is  a 
phenomenon  that  changes  the  ordering  of  natural  modes  with  structural  modification.  The 
phenomenon  often  occurs  between  the  natural  mode  of  interest  and  the  neighboring  orthogonal 
natural  modes.  The  natural  frequency  of  interest  and  the  modal  stiffness  are  drastically  affected, 
resulting  in  convergence  difficulties  and/or  non-optimal  configurations.  In  two-dimensional 
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structures,  Xie  and  his  colleagues  attempted  to  prevent  the  phenomenon  by  increasing  the 
separation  between  frequencies.11  However,  this  method  is  not  applicable  to  three-dimensional 
structures  in  which  the  phenomenon  occurs  unpredictably  during  design  iterations.  In  addition,  as 
the  number  of  elements  eliminated  increases  during  the  iterative  ESO  process,  a  sudden  drop  in 
the  natural  frequency  of  interest  is  observed.  This  is  a  direct  result  of  the  modal  stiffness  decrease 
associated  with  the  current  level  of  element  removal.  The  ESO  method  and  the  associated 
dynamic  sensitivity  number5  fail  to  consider  modal  stiffness  directly.  Consequently,  these 
problems  make  it  harder  for  the  ESO  method  to  detennine  the  Pareto  optimum  in  the 
multi-objective  optimization6  that  contains  dynamic  characteristics. 

In  this  paper,  the  modified  static  control  parameter  using  von  Mises  stress  is  newly 
formulated  to  prevent  non-smooth  changes  in  modal  stiffness  and  the  associated  natural 
frequency  of  interest.  This  modified  static  control  parameter  smooths  the  iterative  process  by 
introducing  correlation  between  the  natural  mode  of  interest  and  its  neighboring  modes,  as  well 
as  a  consideration  of  the  modal  stiffness.  Using  the  proposed  control  parameter,  the  ESO  method 
results  in  a  lightweight  TPS  model  with  improved  static  and  dynamic  characteristics  even  though 
the  mode-switching  phenomena  are  observed  frequently  in  the  iterative  steps.  These 
improvements  are  accomplished  through  efficient  material  volume  reduction. 

In  the  next  section,  control  parameters  for  static  and  dynamic  analysis  based  on  von 
Mises  stress  are  formulated.  Section  III  presents  the  details  of  the  multi-objective  optimization 
technique.  The  evolutionary  optimization  algorithm  is  described  in  Section  IV.  In  section  V,  two 
relevant  metallic  TPS  models  with  different  load  cases  are  presented  to  demonstrate  the 
applicability  of  the  ESO  method  with  the  proposed  control  parameter. 
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3.  SENSITIVITY  ANALYSIS 


Like  most  other  structural  optimization  methods,  the  ESO  method  must  be  iterated  many 
times.  It  is  important  to  search  for  the  removal  positions  where  the  effect  is  high  for  achieving 
the  objective  as  accurately  and  quickly  as  possible.  In  this  section,  two  kinds  of  control 
parameters  (called  sensitivity  numbers5),  based  on  static  analysis  are  presented  as  guidelines  for 
the  removal  of  the  finite  elements.  In  particular,  the  modified  static  control  parameter  is  newly 
developed  to  reduce  the  problems  of  the  weakness  of  the  modal  stiffness  and  the  mode-switching 
phenomenon  in  dynamics. 


3,1  Control  parameter  for  Static  Analysis 

The  stress  level  in  each  element  can  be  measured  by  an  average  of  all  of  the  stress 
components.  For  this  purpose,  the  von  Mises  stress  has  been  one  of  the  most  frequently  used 
criteria  for  isotropic  materials.  The  von  Mises  stress  a™  of  the  /  th  element  for  a 
three-dimensional  structure  is  defined  as 


cr,  = 


~ay.i)2  +(ay,i  +K/  +6(: 


V2 


2  2  2 

T  J  +  T  ,+T 

xy,  l  yz,l  zx 


J  (0 


Here,  crx  ;  ,a  ,,  and  azl  are  normal  stresses  of  the  Ith  element  in  x,  y  ,  and  z  directions, 

respectively,  and  r  ,,  rJZ  , ,  and  z  zx]  are  the  shear  stresses  of  the  /th  element. 

According  to  the  ESO  algorithm,  a  reliable  sign  of  potential  structural  failure  is  excessive 

stress  or  strain,  and  a  low  stressed  element  is  assumed  to  be  under-utilized  in  the  structure. 

Ideally,  the  stress  in  every  part  of  a  structure  should  be  near  the  same  safe  level.  By  gradually 

removing  material  that  has  a  lower  stress,  the  stress  level  in  the  new  designs  become  more  and 

more  unifonn.  By  defining  the  von  Mises  stress  as  the  static  control  parameter  in  the  static 

analysis,  a  number  of  elements  with  the  smallest  static  control  parameter  are  removed  so  that  the 

increase  of  the  maximum  von  Mises  stress  is  minimal.  The  removal  process  must  be  slightly 

12 

modified  when  the  stresses  are  primarily  thermally  induced. 


3.2  Control  parameter  for  Dynamic  Analysis 

The  dynamic  control  parameter  in  the  ESO  method  5'  11  is  derived  from  the  Rayleigh 
quotient  and  is  based  on  the  natural  mode  of  interest.  This  control  parameter  estimates  only  the 
change  of  the  natural  frequency  of  interest.  No  direct  consideration  is  given  to  the  modal 
stiffness.  Also,  when  mode-switching  is  observed  during  the  optimization  process,  the 
characteristic  of  the  dynamic  control  parameter  varies  drastically  due  to  the  change  in  the  natural 
mode  of  interest.  These  effects  can  make  the  modal  stiffness  very  small  bringing  about  a  drastic 
change  in  the  natural  frequency  of  interest  when  a  large  number  of  elements  are  eliminated  from 
the  structure  through  many  iterative  steps. 
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To  address  the  shortcomings  of  the  traditional  ESO  method,  a  static  control  parameter 
using  von  Mises  stress  is  developed,  which  gives  direct  consideration  to  the  modal  stiffness  as  a 
substitute  for  the  conventional  dynamic  control  parameter  based  on  the  Rayleigh  quotient.  In 
addition,  this  static  control  parameter  is  modified  so  that  consideration  is  given  to  adjacent 
natural  modes,  as  well  as  the  natural  mode  of  interest  thus  preventing  a  rapid  change  in  the 
characteristic  of  the  control  parameter. 

The  derivation  of  the  modified  static  control  parameter  for  dynamics  is  as  follows.  The 
equation  of  motion  for  an  undamped  system  is  described  by 

Mf}+MM={/}  <2> 

where  [m]  and  [/<"]  are  the  global  mass  matrix  and  global  stiffness  matrix,  respectively,  and 
M  and  {/}  are  the  global  nodal  displacement  and  nodal  force  vectors,  respectively.  The 
displacement  in  the  spatial  coordinates  of  Eq.  (2)  can  also  be  expressed  using  the  modal 
coordinates  by 


where  N  is  the  total  number  of  degrees  of  freedom,  and  {®}  and  {f}  represent  natural  modes 
and  modal  displacements,  respectively. 


Using  Eq.  (3)  and  the  characteristics  of  generalized  orthogonality,  the  response  in  the 
spatial  coordinates  is  derived  as 


r=i  -to  mr  +kr 


(4) 


where  co  is  the  natural  frequency,  (®  r  }T  [/V/]  {®  }  =  mr  =  1  , 

i/}  =  Ef 

In  this  method,  the  displacement  of  each  nodal  point  is  computed  by  implementing  the 
concept  of  a  virtual  static  displacement  for  each  mode  shape.  The  virtual  static  displacement  is 
chosen  to  be  the  same  displacement  shape  as  that  of  the  mode  shape.  The  displacement  for  each 
mode  shape  is  derived  by  assuming  a>=  0.  The  absolute  value  of  the  displacement  is  given  by 

K  (5) 


Note  that  \kr  ]  =  diag(kl  ,•••,  kN  ) .  From  Eq.  (5),  \kr]  [®]7  {f}  indicates  a  modal  displacement 
that  is  generated  by  the  external  force  {f} .  The  natural  mode  {®  }  can  be  treated  as  a  response 
(displacement)  by  applying  the  external  force  whose  modal  displacement  is  1  for  the  r th  mode 
and  0  for  other  modes.  When  [f.  }  is  the  external  force  satisfying  this  condition,  {f  }  is  easily 
obtained  by  the  relationship  between  the  displacement  and  external  force,  as  shown  in  Eq.  (6): 
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(6) 

In  general,  the  magnitude  of  the  external  force  will  vary  from  mode  to  mode,  with  higher 
natural  modes  requiring  larger  external  forces  due  to  the  complexity  of  the  higher  mode  shapes. 
This  tendency  is  exploited  in  the  development  of  a  novel  control  parameter  which 
simultaneously  considers  all  frequencies  of  interest.  The  hypothesis  uses  relative  force 
magnitudes  as  a  surrogate  for  “closeness”  of  natural  frequencies.  The  nodal  force  vectors  are 
scaled  such  that  large  nodal  forces  (corresponding  to  complicated  mode  shapes)  are  reduced  by 

the  factor  (]|^||/||Fr||)P  .  Thus  {Fr}new  is  given  by  Eq.  (7): 


K 


IK 

IK 


K 


(7) 


Here,  subscript  i  shows  the  order  of  interest,  ( 0  <  i  «  N  ).  Exponent  P,  which  adjusts  the 
magnitude  of  the  scaling  factor  between  the  i th  natural  mode  and  the  r th  natural  mode,  is  the 
correlation  factor  (typically  P  =  1).  A  larger  value  of  P  emphasizes  the  relative  importance  of  the 
natural  mode  of  interest.  The  same  scaling  can  be  used  for  the  response 
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The  scaled  eigenvectors  {® r  }new  are  then  applied  as  static  displacement  fields. 

Using  Eqs.  (1)  and  (8),  von  Mises  stresses  of  each  element  are  obtained  from  the  z th 
natural  mode  to  the  n  th  natural  mode.  Each  natural  mode  is  used  as  a  displacement  vector,  and 
the  stress  components  for  each  natural  mode  are  obtained  by  using  the  strain-displacement  and 
stress-strain  relations.  Finally,  for  each  finite  element,  the  maximum  von  Mises  stress  value 
between  the  i th  and  the  n  th  natural  modes  is  selected  as  the  new  dynamic  control  parameter,  as 
expressed  by  Eq.  (9). 


crv,m  =max 


cr;,  /  >  °’i+ 1,  I  I&i+l,  l  >' '  ■>CrH, 


(0 


<  i  «  n,  n  < 


N) 


(9) 


Since  Eq.  (9)  considers  not  only  the  z'th  natural  frequency,  but  also  its  neighboring  natural 
frequencies,  the  smooth  change  of  natural  frequency  can  be  possible  even  if  mode-switching 
occurs  in  the  iterative  processes.  For  example,  when  other  natural  frequencies  are  much  higher 
than  the  i th  natural  frequency,  only  the  maximum  von  Mises  stress  of  the  i th  natural  mode  is 

selected  for  each  finite  element  because  (j| Fi  ||/|| Fr  ||)P  for  other  natural  modes  has  small  values 

(«1).  On  the  other  hand,  as  the  neighboring  natural  frequencies  approach  that  of  the  zth  natural 
frequency  during  the  iterative  process,  other  natural  modes  are  also  gradually  considered 

(i.e.,(j|if  |/||F(.  ||)P  =  1).  Note  that  it  is  not  necessary  to  consider  lower  frequencies  than  the  z'th 
natural  frequency  because  they  don’t  switch  with  the  i th  natural  frequency  in  most  cases. 
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Since  the  proposed  dynamic  control  parameter  is  derived  by  using  the  von  Mises  stress,  a 
number  of  elements  having  some  of  the  smallest  von  Mises  stresses,  that  is,  the  most  inefficient 
elements  are  removed  from  the  structure,  identical  to  the  static  optimization  problem  of  the  ESO 
method. 
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4.  MULTI-OBJECTIVE  OPTIMIZATION  TECHNIQUE 


In  the  weighted  objectives  method,  this  multi-objective  problem  is  converted  to  a 
single-objective  problem  by  using  a  weighted  sum  of  the  original  multiple  objectives.  That  is,  the 
control  parameter  for  the  static  characteristic  in  Eq.  (1)  and  the  control  parameter  for  the 
dynamic  characteristic  in  Eq.  (9)  are  combined  together  to  form  a  new  single  criterion: 

Fr“=W,Ru+W2R2:,='ZwjRi,,  (10) 

7=1 

where  F"'"1 "  is  the  multi-objective  function  that  detennines  element  removal  for  each  element  /. 

2 

Wj  is  the  j th  criterion  weighting  factor  with  0  <  IT  <  1  and  ^  IT  =  1 .  R  j ,  =  <j ™l  j <j J"max  is 

7=1 

the  ratio  of  the  j th  criterion  control  parameter  o  v"\  for  the  Ith  element  to  the  maximum  value 
of  the  j th  criterion  control  parameter  cr/y'"max .  Since  both  control  parameters  are  derived  from 
von  Mises  stress,  a  number  of  elements  with  the  smallest  F'’rul“  are  removed  so  that  both 
characteristics  can  be  improved,  as  in  the  conventional  ESO  method  for  statics. 
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5.  EVOLUTIONARY  STRUCTURAL  OPTIMIZATION 

ALGORITHM 

The  ESO  algorithm  of  the  weighted  objectives  method  for  stress  and  the  fundamental 
natural  frequency  is  implemented  as  follows: 

i)  Discretize  the  structure  using  a  fine  mesh  of  finite  elements. 

ii)  Solve  the  linear  static  analysis  problem  under  thermal  loads  (or  mechanical  loads). 

iii)  Calculate  the  von  Mises  stress  by  using  Eq.  ( 1 ). 

iv)  Solve  the  eigenvalue  problem. 

v)  Calculate  the  von  Mises  stress  by  using  Eq.  (9). 

vi)  Combine  the  two  control  parameters  by  using  Eq.  (10). 

vii)  Remove  several  elements  from  the  current  structure  that  have  relatively  small 
contributions. 

viii)  Repeat  steps  ii)  to  vii)  until  the  weight  of  the  structure  reaches  the  predetermined  value. 

The  typical  ESO  process  is  operated  by  a  rejection  ratio  (RR)  that  determines  which 
elements  to  eliminate.  During  the  iterative  process,  the  rejection  ratio  is  gradually  increased  by 
the  evolutionary  rate  ( ER  ).5  The  rejection  ratio  is  simply  implemented  by  comparing  the  relative 
reference  levels  of  candidate  elements  with  the  maximum  value.  In  general,  the  maximum  value 
of  the  reference  criterion  doesn't  decrease  during  the  evolutionary  process  for  most  mechanical 
load  cases.  This  allows  the  threshold  levels  to  steadily  increase  with  material  removal. 
Conversely,  for  cases  involving  thennal  stress,  the  maximum  stress  level  can  decrease  as  the 
evolutionary  process  progresses.  "  For  this  reason,  a  new  rejection  ratio  RR„ew  is  determined  as 

the  a  th  lowest  F muIn  value.  Note  that  ex  should  be  sufficiently  small  to  guarantee  a  smooth 
change  between  two  iterations.  In  this  method,  a  is  set  to  1~2%  of  the  total  number  of 
elements  for  the  initial  iteration. 
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6.  THERMAL  PROTECTION  SYSTEMS  DESIGN 


Because  the  TPS  comprises  the  external  surface  of  a  spacecraft,  several  requirements 
must  be  satisfied  to  design  a  TPS  for  all  environments  experienced  by  the  vehicle: 

a)  In  order  to  minimize  operational  costs,  a  TPS  should  be  as  lightweight  as  possible 
while  maintaining  its  stiffness  and  high  natural  frequency. 

b)  A  thin  TPS  plate  is  needed  to  protect  a  spacecraft  from  heat,  and  a  frame  structure 
may  be  attached  to  prevent  fluttering  due  to  aerodynamic  and/or  acoustic  loading.  If 
required,  a  honeycomb  sandwich  can  be  affixed  to  the  plate  if  the  spacecraft 
encounters  a  maximum  surface  temperature  for  an  extended  period  of  time. 

c)  A  long  support  connecting  the  plate  to  the  fuselage  is  recommended  to  decrease  heat 
transfer  to  the  fuselage.  This  construction  provides  for  extra  volume  of  insulation  to 
be  placed  in  the  cavity,  and  facilitates  the  generation  of  novel  attachment  by  ESO  (or 
other  topology)  methods. 

d)  The  maximum  thennal  stress  in  the  TPS  is  required  to  remain  below  yield  strength  to 
avoid  permanent  defonnation. 

In  this  section,  two  relevant  TPS  models  with  different  initial  conditions  are  presented  to 
demonstrate  the  effectiveness  of  the  ESO  method  with  the  proposed  control  parameter.  An 
efficient  means  for  designing  a  TPS  is  shown  for  various  requirements  by  employing  the 
proposed  control  parameter.  The  initial  TPS  design  is  shown  in  Fig.  1.  Inconel  alloy  693  is 
utilized  because  of  its  excellent  resistance  to  metal  dusting  and  high  temperature  corrosion,  as 
well  as  for  its  favorable  fabrication  and  joining  properties.  The  TPS  models  are  discretized  with 
hexahedron  isoparametric  linear  elements  for  finite  element  analysis,  and  the  p-version  of  the 
linear  finite  element  method  with  Guyan  reduction  is  applied  to  enhance  the  accuracy  of  the 
analysis  and  reduce  computational  effort.14  Structural  characteristic  matrices  with  more  degrees 
of  freedom  are  generated,  although  the  sizes  of  the  matrices  are  identical  to  those  obtained  with 
the  h-version  finite  element  method.  For  both  examples,  a  ,  a  parameter  for  the  element  removal, 
as  shown  in  the  previous  section,  is  set  at  1%,  and  n  in  Eq.  (9)  is  set  at  10. 


6.1  Example  1:  TPS  Model  with  External  Forces 

The  conventional  ESO  method  and  the  new  ESO  method  with  the  proposed  control 
parameter  are  both  applied  to  a  TPS  model  with  a  mechanical  load  case.  The  TPS  model,  which 
has  the  dimensions  of  A=0.5  m,  B=0.45  m,  C=0.03  m,  and  D=0.27  m,  as  shown  in  Fig.  1,  is 
loaded  at  room  temperature,  with  a  Young's  modulus  of  196  GPa,  a  material  density  of  7770 
kg/nr,  and  a  Poisson's  ratio  of  0.32.  The  model  is  meshed  with  0.05  mx0.05  mx0.03  m 
rectangular  isoparametric  elements  with  fixity  prescribed  on  the  bottom  side.  External  tractions 
of  magnitude  10  N  in  the  X  direction  are  applied  to  each  nodal  point  on  the  top  side.  The 
elements  identified  with  the  external  force,  which  include  those  in  the  plate-frame  region,  are  set 
as  unremovable  elements.  This  unremovable  region  renders  the  evolutionary  process  unstable  in 
the  conventional  ESO  method  as  adjacent  regions  with  high  modal  mass  are  eliminated  without 
consideration  of  the  modal  stiffness. 

The  optimization  problem  is  to  minimize  the  maximum  von  Mises  stress  and  to  maximize 
the  fundamental  natural  frequency  while  reducing  the  TPS  weight.  These  two  objectives  are 
applied  as  weighted  objectives  as  shown  in  Eq.  (10).  The  relationship  between  maximum  von 
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Mises  stress  and  the  fundamental  natural  frequency  is  investigated  by  varying  both  weighting 
factors  and  volume  reduction. 

Figure  2  shows  the  comparison  between  the  fundamental  natural  frequency  and  the 
maximum  von  Mises  stress  when  reductions  in  volume  are  20%  and  50%.  The  plots  show  results 
for  both  the  conventional  ESO6  and  the  proposed  method.  The  weighting  factor  for  the  dynamic 
characteristic  (Wd  ic)  and  the  weighting  factor  for  the  static  characteristic  (Wstatic)  for  each 

point  in  Fig.  2  are  as  follows: 


(A)  and  (A’): 

W 

dynamic 

■  w 

*  static 

=  0. 0:1.0, 

(B)  and  (B’): 

w 

dynamic 

w  . 

static 

=  0.2:0. 8, 

(C)  and  (C): 

w 

dynamic 

w  . 

static 

=  0.5:0. 5, 

(D)  and  (D’): 

w 

dynamic 

w  . 

static 

=  0. 8:0.2, 

(E)  and  (F): 

w 

dynamic 

w  . 

static 

=  1. 0:0.0. 

Points  of  (A)-(E)  and  (A')-(E')  are  results  from  the  conventional  and  the  proposed  method, 
respectively. 

In  the  proposed  method,  it  is  observed  that  each  objective  is  improved  by  increasing  its 
respective  weighting  factor  (e.g.  both  the  fundamental  natural  frequency  and  the  maximum  von 
Mises  stress  increase  when  the  dynamic  weighting  factor  is  increased  and  the  static  weighting 
factor  is  decreased).  Because  the  proposed  method  applied  to  the  modified  static  control 
parameter  directly  considers  the  mode-switching  phenomenon  and  the  modal  stiffness,  smooth 
changes  are  observed  in  the  static  and  dynamic  characteristics  even  as  a  large  number  of 
elements  are  removed  through  many  iterations.  It  is  shown  that  any  improvement  of  one  criterion 
requires  a  clear  tradeoff  with  the  other,  and  a  clear  Pareto  solution  (or  a  Pareto  curve)  can  be 
observed. 

Unlike  the  proposed  method,  a  desirable  solution  is  not  obtained  in  the  conventional 
method.  As  shown  in  Fig.  2,  the  maximum  von  Mises  stresses  obtained  with  the  conventional 
method  are  higher  than  those  obtained  from  the  proposed  method  with  increasing  Wdynamic .  Even 

though  the  conventional  method  addresses  the  relative  ratio  between  modal  stiffness  and  modal 
mass,  the  deficiency  arises  from  no  direct  consideration  of  modal  stiffness.  In  addition,  it  is 
shown  that  the  fundamental  natural  frequencies  at  Wdynamic :  Wstatic=\  .0:0.0  are  decreased  by  this 

lack  of  modal  stiffness  as  well  as  the  changing  control  parameter  characteristics  due  to 
mode-switching.  Determining  the  trend  of  the  Pareto  solution  with  volume  reduction  is 
extremely  difficult.  Consequently,  compared  with  the  conventional  method,  the  proposed  control 
parameter  for  an  eigenvalue  problem  can  be  used  to  design  an  optimum  structure  that  requires 
the  elimination  of  a  large  number  of  elements  through  many  iterative  processes  in  addition  to 
favorable  dynamic  modification. 

6.2  Example  2:  TPS  Model  with  Heat  Transfer  Considerations 

When  a  spacecraft  re-enters  the  atmosphere  from  Low  Earth  Orbit,  a  metallic  TPS  creates 
a  temperature  distribution  due  to  tremendous  heat  flux.  The  temperature  distribution  changes  the 
material  characteristics  of  a  metallic  TPS,  such  as  thennal  conductivity,  Young's  modulus,  etc. 
This  temperature  dependence  in  the  thermal  problem  results  in  nonlinearity  and  requires  an 
iterative  solution  technique. 
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For  the  second  example,  a  transient  conduction  analysis  is  applied  to  the  initial  TPS 
model  in  Fig.  1.  It  is  assumed  that  the  [void]  space  inside  the  TPS  and  the  vertical  surface 
boundaries  (with  the  exception  of  the  surface  of  the  plate)  are  insulated.  The  plate  temperature 
(Tpiate)  on  the  top  side  varies  according  to  a  temperature  profde,  as  shown  in  Ref.  15.  In  this 
method,  the  temperature  is  simply  represented  as  Tplate  =  1 .2540t  +  200  (°C)  when  t  <  500 

(sec.),  and  T  late  =  827  (°C)  when  t  >  500  (sec.),  regardless  of  the  elimination  of  elements 

from  the  structure.  This  assumption  is  a  severe  restriction.  In  the  true  structure,  convection  and 
radiation  play  major  roles  in  the  transfer  of  energy.  For  instance,  to  address  the  internal  radiative 
heat  transfer  that  will  occur  as  material  is  removed,  an  efficient  method  for  constructing  view 
factors  during  the  iterative  process  would  be  needed.  Since  the  purpose  of  this  work  was  to 
develop  an  ESO  method  for  simultaneous  consideration  of  thennal  stress  and  frequency, 
prescribed  temperature  boundary  conditions  were  used  to  ensure  tractability. 

A  numerical  time  integration  scheme,16  similar  to  Newmark's  method,  is  used  to  solve  the 
following  transient  heat  transfer  problem 


[M,  ]  +  P  [k, ])  {tm }  =  [" T  [M,  ]- (1  -  P)[K, } 


At 


At 


P){F,}+ p{fm 


(ii) 


where  \M t  ] :  Consistent-mass  matrix 

\Kt  ] :  Thennal  conductivity  matrix 

{if  }:  Force  matrix  defined  by  heat  source  (or  heat  sink),  heat  flux,  and  convection 

At :  Time  step 

(3  :  Parameter,  which  is  set  as  0.5  by  the  Crank-Nicolson  rule 

{T. } ,  {Tm  } :  Temperatures  at  time  f  and  tM  ,  respectively 


To  reduce  the  computational  time  of  the  transient  analysis,  only  the  tendency  of  the  temperature 
profile  at  t =1500  (sec.)  is  evaluated  by  assuming  that  [mJ  and  [/<",  ]  as  constant  with  time  as 

thennal  conductivities  ( Kxx  =  Kyy=Kzz=  20  W/(m-°C)),  density  ( p  =1110  kg/m3),  and  specific 

heat  (c=530  J/(kg-°C))  are  treated  as  temperature-independent  variables.  A  constant  time  of 
At  =  10  sec.  is  considered.  The  initial  temperature  distribution  {T()  }  for  each  iterative  process  is 

assumed  to  be  a  linear  profile  through  the  thickness  with  127  °C  on  the  top  side  and  27  °C  on  the 
bottom  side. 


6.2.1  Availability  of  the  Proposed  Control  Parameter  in  the  Heat  Transfer  Problem 

Both  conventional  and  proposed  methods  are  compared  to  demonstrate  the  applicability 
of  the  proposed  control  parameter  for  topology  optimization  of  a  TPS.  The  TPS  model  has  the 
dimension  of  A=0.48  m,  B=0.42  m,  C=0.03m,  and  D  is  specifically  set  as  0.27  m  longer  than 
conventional  structures17  to  reduce  thermal  stress  in  the  fixed  regions  at  the  bottom  side 
(interface  to  low-temperature  material).  The  plate  and  frame  regions  are  treated  as  unremovable. 
The  model  is  divided  into  0.03  mx0.03  mx0.03  m  regular  hexahedron  isoparametric  elements, 
and  the  0.06  mx0.06  m  bottom  corner  regions  are  considered  fixed.  The  coefficient  of  expansion 
and  Young's  modulus  are  considered  as  functions  of  temperature.13  In  much  the  same  way  as 
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Example  1,  the  weighted  objectives  method  considers  both  the  maximum  thennal  stress  and  the 
fundamental  natural  frequency  by  changing  both  weighting  factors  and  volume  reduction. 

Figure  3  shows  the  changes  in  fundamental  natural  frequencies  and  maximum  thermal 
stresses  by  varying  weighting  factors  and  sensitivities.  Both  methods  improve  the  fundamental 
natural  frequencies  by  increasing  the  weighting  factors  for  dynamics,  as  in  the  mechanical  load 
case.  In  general,  when  a  structure  is  over-designed,  in  other  words,  when  a  structure  consists  of 
more  elements  than  needed,  the  fundamental  natural  frequency  can  be  increased  by  removing 
inefficient  elements.  However,  when  the  stiffness  of  the  structure  becomes  weak,  the 
fundamental  natural  frequency  is  decreased.  In  Fig.  3,  as  a  large  number  of  elements  are  removed 
from  the  structure  through  many  iterative  steps,  the  proposed  method  keeps  the  fundamental 
natural  frequencies  much  higher  than  that  obtained  from  the  conventional  method,  even  though 
there  is  no  inertia  consideration  in  the  proposed  control  parameter.  Conversely,  the  results  using 
the  conventional  control  parameter  show  that  the  fundamental  natural  frequency  decreases 
quickly  due  to  (1)  the  small  modal  stiffness  in  the  connecting  region  between  the  assumed  fixity 
and  the  plate-frame  region  and  (2)  the  characteristic  change  in  the  control  parameter  due  to  the 
discontinuous  change  in  the  natural  modes  (i.e.,  mode-switching  phenomenon). 

On  the  other  hand,  the  analysis  of  the  change  in  the  maximum  thermal  stress  reveals 
different  behavior  when  compared  with  the  mechanical  load  case.  A  higher  static  weighting 
factor  does  indeed  result  in  a  lower  maximum  thermal  stress.  In  the  conventional  method,  the 
minimization  of  the  maximum  thennal  stress  can  be  obtained  by  applying  a  large  weighting 
factor  for  the  static  control  parameter  in  the  early  stages  of  the  analysis,  similar  to  the  mechanical 
load  case.  As  the  volume  removal  percentage  increases,  the  maximum  thennal  stress  that  results 
from  dynamic  considerations  alone  becomes  less  than  the  maximum  thennal  stress  that  results 
from  static  considerations  alone.  These  seemingly  inconsistent  results  are  due  to  heat  transfer 
considerations.  With  constant  temperature  boundary  conditions  prescribed,  the  heat  flux  to  the 
bottom  side  is  decreased  with  the  elimination  of  elements  that  connect  the  plate-frame  to  the 
support  region.  The  results  from  the  proposed  method  have  analogous  tendencies  to  the 
conventional  method  with  changing  weighting  factors.  However,  the  proposed  method  can 
restrain  the  maximum  stresses  for  Wdynamic  :  Wstatic  =1.0: 0.0  from  increasing  at  early  stages  of 

the  process. 

Figure  4  shows  resulting  models  close  to  900  Hz  in  fundamental  natural  frequency  by 
both  the  conventional  and  proposed  methods.  Each  model  is  developed  using  solely  dynamic 
considerations,  that  is,  Wd  mic  :  Wstatic  =  1.0  : 0.0  .  The  fundamental  natural  frequency,  maximum 

stress,  and  the  number  of  elements  result  in  933.6  Hz,  0.523  GPa,  and  740  (155.2  kg)  using  the 
conventional  method,  and  906.8  Hz,  0.330  GPa,  and  568  (119.2  kg)  by  the  proposed  method.  The 
main  difference  between  the  two  resulting  topologies  is  the  location  of  the  elements  that  are 
adjacent  to  the  plate-frame  region:  one  is  located  towards  the  center  of  the  plate-frame  region, 
and  the  other  is  located  towards  the  edges  of  the  plate-frame  region. 

From  the  above  results,  the  ESO  method  with  the  proposed  control  parameter  can  be  used 
to  design  a  lightweight  TPS  model  with  a  high  fundamental  natural  frequency  and  a  low 
maximum  thermal  stress.  Using  the  assumption  that  only  conduction  affects  the  support  region,  a 
practical  TPS  design  is  developed  by  using  Wdynamic  :  Wstatic  =1.0: 0.0  in  the  next  section. 
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6.2.2  TPS  Design  with  Heat  Transfer  Problem 

In  contrast  to  the  static  case,  the  convergence  of  the  eigenvalue  maximization  problem 
presents  a  more  daunting  task.  In  an  evolutionary  optimization  that  requires  a  large  number  of 
element  removals  through  many  iterative  steps,  the  ordering  of  natural  modes  can  be  altered  by 
the  selection  of  the  initial  geometry.  This  input  can  result  in  considerable  variations  to  the 
problem’s  structural  characteristics  (i.e.  neighboring  natural  frequencies,  weight  of  the  structure, 
etc.).  Therefore,  based  on  the  topology  results  of  the  previous  section  (Fig.  4  (b)),  the  initial  TPS 
model  is  chosen  to  be  that  of  Fig.  4  (b)  with  a  slight  modification.  Element  are  removed  from 
inside  the  structure,  as  shown  in  Fig.  5.  The  shape  optimization  method  using  the  ESO  technique 
(called  ’nibbling  ESO’5)  is  then  applied  to  the  support  region  to  avoid  the  checkerboard  pattern 
that  occurred  in  Fig.  4  (b).  That  is,  only  the  structural  boundaries  in  the  Z  direction  between 
Layer  1  and  Layer  10  (with  the  exception  of  the  unremovable  region  of  Layer  0)  are  removed 
from  the  model.  The  frame  region,  as  well  as  the  support  region,  is  designed  to  reduce  the  TPS 
weight.  Using  the  proposed  method,  the  objectives  take  the  following  form: 

Minimize  TPS  weight 
subject  to 

1 .  maximization  of  the  fundamental  natural  frequency  (  fx  >  900  Hz) 

2.  minimization  of  the  maximum  thermal  stress  below  yield  strength 
(amax<0.3GPa) 

In  early  iterations,  the  fundamental  natural  mode  is  the  bending  mode  of  the  support 
region.  The  elements  of  the  frame  part  are  eliminated  to  improve  the  dynamic  characteristics  of 
the  support.  The  local  bending  mode  of  the  thin  top  plate  becomes  the  lowest  mode  of  the 
structure.  If  the  original  model  was  too  coarse  or  too  much  material  was  removed  at  each 
iteration,  the  plate  may  not  satisfy  the  frequency  constraint  and  could  result  in  undesirable  flutter 
or  acoustic  excitation.  The  problem  becomes  a  multi-scale  phenomenon.  Inefficient  material  still 
exists  in  the  base  but  the  optimization  algorithm  is  focusing  on  the  low  frequency  of  the  plate. 
Therefore,  to  address  the  problem  efficiently,  the  design  optimizations  for  the  support  region  and 
the  frame  region  are  conducted  separately. 

TPS  Support  Design 

First,  the  support  region  for  the  initial  model  in  Fig.  5  is  set  up  with  identical  boundary 
conditions  and  mesh  size,  as  with  that  of  previous  section  V.B.l.  However,  additional  fixed 
regions  are  established  at  the  edges  of  the  plate  to  prevent  a  local  mode  from  occurring  in  the 
unremovable  region  due  to  dynamic  considerations.  By  including  additional  stiffness  [kx,  kY,  kz] 
=  [0,  0,  10']  (N/m),  the  local  mode  of  the  unremovable  region  will  have  a  higher  natural 
frequency  when  compared  with  the  fundamental  natural  frequency  of  the  support's  first  bending 
mode.  Figure  6  shows  the  change  in  the  fundamental  natural  frequency  and  the  maximum 
thennal  stress  as  a  function  of  volume  reduction,  which  is  the  ratio  percentage  between  the 
current  and  the  fully-populated  structure's  volume.  For  this  case,  the  structure  at  77.7%  volume 
removal  is  the  most  lightweight  that  satisfies  the  two  constraints  (Fig.  7).  When  additional 
stiffness  is  removed  from  the  plate,  the  fundamental  natural  frequency  is  less  than  900  Hz  due  to 
the  increase  in  modal  mass  of  the  unremovable  region.  In  this  case,  the  TPS  model  has  the 
fundamental  natural  frequency  of  871.9  Hz,  and  the  maximum  thermal  stress  occurs  in  the 
support  region. 
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TPS  Frame  Design 

Using  the  modified  TPS  model  in  Fig.  7,  the  model  is  then  designed  to  be  lightweight 
with  the  fundamental  natural  frequency  greater  than  900  Hz  by  modifying  the  frame  shape.  In  the 
frame  design,  it  is  possible  to  increase  the  fundamental  natural  frequency  because  all  elements  in 
the  frame  region  have  a  high  fundamental  modal  mass.  The  modal  stiffness  of  the  plate  and  the 
frame  can  be  also  considered  by  applying  the  proposed  control  parameter.  The  plate  and  the 
frame  regions  are  remeshed  with  two  element  sizes,  0.03  mx0.03  mx  0.007  m  in  the  frame  and 
0.03  mx0.03  mx  0.002  m  in  the  plate.  As  in  the  element  size  of  the  plate,  the  thickness  is  set  at 
0.002  m.  The  four  corners  at  the  bottom  of  the  support  region  are  assumed  fixed  as  in  the 
example  of  Section  V.B.l.  Shape  optimization  is  applied  by  eliminating  elements  from  the 
bottom  surface  of  the  frame  region  for  ease  of  manufacture,  as  with  the  support  design.  Because 
the  maximum  thermal  stress  exists  in  the  support  region,  no  thermal  stress  analysis  is  applied. 
Only  the  temperature  profile  of  the  TPS  model  in  Fig.  7  is  utilized  for  the  evaluation  of  the 
material  properties. 

Figure  8  shows  the  evolutionary  history  of  the  fundamental  natural  frequency  due  to 
modifications  of  the  frame  region.  Up  to  the  41st  iteration,  the  fundamental  frequency  increases 
because  excessive  modal  mass  is  removed  from  the  frame.  Then,  the  frequency  drastically 
decreases  because  the  modal  stiffness  is  reduced  by  the  lack  of  elements  in  the  frame  region, 
even  though  the  frame  has  excessive  modal  mass.  In  the  frame  design,  the  modified  TPS  model 
satisfies  the  frequency  constraints  at  the  45  th  iteration,  and  has  a  fundamental  natural  frequency 
of  920.2  Hz. 

Finally,  a  transient  heat  transfer  analysis  is  applied  to  the  TPS  model  at  the  45th  iteration 
for  the  evaluation  of  the  exact  solution.  In  this  case,  the  model  has  the  fundamental  natural 
frequency  of  919.8  Hz,  the  maximum  thermal  stress  of  0.228  GPa,  and  816  elements  (76.50  kg), 
which  consist  of  224  elements  (3.13  kg)  in  the  plate,  316  elements  (15.47  kg)  in  the  frame,  and 
276  elements  (57.90  kg)  in  the  support  region.  The  fundamental  natural  frequency  is  reduced 
slightly  by  the  change  of  the  material  properties,  and  the  maximum  thennal  stress  is  also  reduced 
slightly  by  the  volume  reduction  of  the  frame.  The  final  modified  TPS  model  and  the 
cross-section  of  the  frame-plate  are  shown  in  Fig.  9.  To  prevent  flutter  of  the  plate,  it  is  shown 
that  the  stiffness  of  the  center  region  on  the  plate  is  reinforced  by  the  existence  of  the  elements  in 
the  frame,  even  though  the  elements  have  a  high  modal  mass  about  the  natural  mode  of  the  plate. 
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7.  CONCLUSIONS 


A  metallic  TPS  concept  provides  the  possibility  of  decoupling  the  thennal  and  structural 
functions  by  providing  metallic  substructures  to  encapsulate  internal  insulation,  maintain  panel 
shape,  and  support  thennal  and  mechanical  loads.  This  decoupling  allows  the  use  of  structurally 
efficient  internal  insulations.  However,  these  two  functions  can  not  be  decoupled  perfectly,  so 
each  substructure  in  the  metallic  TPS  model  should  be  designed  to  reduce  thermal  defects,  as 
well  as  to  allow  for  structural  modification. 

In  this  paper,  a  multi-objective  optimization  problem  for  the  thennal  stress  and  the 
fundamental  natural  frequency  was  conducted  to  make  a  lightweight  TPS  model  by  using  the 
Evolutionary  Structural  Optimization  (ESO)  algorithm.  Two  objectives  were  optimized 
simultaneously,  namely,  the  maximization  of  the  fundamental  natural  frequency  and  the 
minimization  of  the  maximum  thermal  stress,  through  efficient  volume  reduction.  Specifically, 
the  modified  static  control  parameter  based  on  static  analysis  was  newly  proposed  for  the  TPS 
design  concerned  with  dynamic  analysis.  Through  the  comparison  between  the  conventional  and 
the  proposed  control  parameter,  it  was  shown  that  the  proposed  control  parameter  can  prevent 
several  problems,  which  are  induced  by  failing  to  consider  the  mode-switching  phenomenon  and 
the  reduced  modal  stiffness  directly  when  a  large  number  of  elements  are  eliminated  through 
many  iterative  steps  (although  the  proposed  control  parameter  ignores  inertia  terms  of  the 
equation  of  motion).  Additionally,  an  efficient  way  to  obtain  a  metallic  TPS  model  was  shown  by 
designing  the  frame  and  the  support  region  separately  using  the  ESO  method  with  the  proposed 
control  parameter.  In  future  work,  consideration  will  be  given  to  more  realistic  thennal 
conditions.  Inclusion  of  convection  and  radiation  as  well  as  a  maximum  temperature  constraint 
will  potentially  result  in  new  and  innovative  designs. 
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Fig.  2  Relationship  between  the  fundamental  natural  frequency  and  the  maximum  stress 
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(a)  Frequency  histories  by  the  conventional  method 


(c)  Max.  stress  histories  by  the  conventional  method 


(b)  Frequency  histories  by  the  proposed  method 


(d)  Max.  stress  histories  by  the  proposed  method 


Fig.  3  Evolutionary  histories  of  the  fundamental  natural  frequencies  and  the  maximum 

thermal  stress 
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(a)  Modified  TPS  model  by  the  conventional  method 
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(b)  Modified  TPS  model  by  the  proposed  method 


Fig.  4  Resultant  TPS  models  with  900  Hz  fundamental  natural  frequency 


21 


Z  <m) 


LAYER  0 
LAYER  1 
LAYER 2 
LAYER 3 
LAYER 4 


LAYER 9 
LAYER  10 


Fig.  5  Initial  model  for  designing  TPS  support 
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Fig.  6  Evolutionary  histories  for  TPS  support  design 
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Fig.  7  Optimum  TPS  support 
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Fig.  8  Evolutionary  history  of  the  fundamental  natural  frequency  for  TPS  frame  design 
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(a)  Optimum  TPS  model 
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(b)  Cross-section  in  the  plate-frame  region 


Fig.  9  Optimum  TPS  model  including  heat  transfer  effects 
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=  global  nodal  force  vector 
=  external  force  vector  for  the  r  th  natural  mode 
=  global  stiffness  matrix 
=  global  mass  matrix 
=  number  of  degrees-of-freedom 
=  total  number  of  degrees-of-freedom 
=  correlation  factor 
=  ratio  between  cr  ™  and  cr™max 

=  multi-objective  control  parameter  for  the  I th  element 
=j th  criterion  weighting  factor 
=  global  nodal  displacement  vector 
=  global  nodal  acceleration  vector 
=j th  von  Mises  stress  for  the  I th  element 
=  maximum  value  of  the  j th  von  Mises  stress 

=  parameter  for  the  element  removal 
=  natural  frequency 
=  natural  mode 
=  modal  displacement  vector 
=  L2  nonn 
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